function [rstring, xlabel, xlabelc] = tdlocal(rtype, xc, yc, zc, cs)
% tdlocal  - lookup local database file for TalLabel
%
%       FORMAT:     tdlabel = tdlocal(rtype, x, y, z, size)
%  or
%       FORMAT:     tdvoi   = tdlocal(rtype, voiname)
%
% Input fields:
%       rtype       numeric type selection: 2: TalLabel, 3: TalCube, 5: NGM
%       x,y,z       Talairach coordinates of request
%       size        only required for types 3 and 5; where with type 3
%                   this is the width, height, and depth of the searched
%                   cube, while with type 5 its a maximum distance;
%                   for type 3, its maximum is 13, for type 5 it's 8.5
%
%       voiname     string that is matched against the list of labels
%
% as with the original network database, tdlocal returns its reply in one
% continuous row. hence, it's preferable to use tdclient.m instead for
% requests.
%
% for an extension, NGM search does **not** search in a cubed pattern, but
% rather spherically, hence it is far more accurate than the same function
% in other localized TD client versions.
%
% See also tdclient.

% Version:  v0.5c
% Build:    6120415
% Date:     Dec-04 2006, 3:15 PM CET
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% enough arguments ?
if nargin < 2 || ...
   ~isa(rtype, 'double') || ...
    numel(rtype) ~= 1 || ...
    isinf(rtype) || ...
    isnan(rtype) || ...
    rtype < -6 || ...
    rtype > 6 || ...
   (~ischar(xc) && nargin < 4)
    error( ...
        'BVQXtools:TooFewArguments',...
        'Bad or too few arguments.' ...
    );
end

% persistent variable for internal database representation
persistent tbase_opts;
if isempty(tbase_opts)

    % db is not loaded yet, loading (and saving) should be tried though
    tbase_opts.dbloaded = false;
    tbase_opts.reload   = true;

    % database filenames
    tbase_opts.mfile    = 'database.mat';
    tbase_opts.dfile    = 'database.dat';
    tbase_opts.lfile    = 'database.txt';

    % x,y,z size
    tbase_opts.xsize    = 144;
    tbase_opts.xtrans   = 73;
    tbase_opts.ysize    = 176;
    tbase_opts.ytrans   = 105;
    tbase_opts.zvals    = [-40:4:-4, -1, 1, 4:4:32, 35:5:65];
    tbase_opts.zsize    = length(tbase_opts.zvals);

    % populate zlookup
    zlookupb = tbase_opts.zvals;
    zlookupi = 1:length(zlookupb);
    tbase_opts.zspat    = [0 1 -1 2 -2];
    for lc = 2:length(tbase_opts.zspat)
        ztest = tbase_opts.zvals - tbase_opts.zspat(lc);
        for ic = 1:length(ztest)
            if ~any(zlookupb == ztest(ic))
                zlookupb(end+1) = ztest(ic);
                zlookupi(end+1) = ic;
            end
        end
    end
    [zlookupb, zlsi] = sort(zlookupb);
    zlookupi = zlookupi(zlsi);
    tbase_opts.zlub = zlookupb;
    tbase_opts.zlui = zlookupi;
    tbase_opts.zmax = max(tbase_opts.zlub);
    tbase_opts.zmin = min(tbase_opts.zlub);

    % full path and mat filename
    tdsrcfld = bvqxtools_path('tal');
    tdmatfile = [tdsrcfld '/' tbase_opts.mfile];

    % preloading file exists and should be used
    if exist(tdmatfile, 'file') == 2 && ...
        tbase_opts.reload

        % try loading
        try
            saved = struct;
            load(tdmatfile);

            % and get contents
            tbase_opts.labels    = saved.labels;
            tbase_opts.lookup    = saved.lookup;
            tbase_opts.dbloaded  = true;
            clear saved;
        catch
            % nothing
        end
    end

    % db not yet loaded
    if ~tbase_opts.dbloaded

        % get full filenames
        tdsrcfile = [tdsrcfld '/' tbase_opts.dfile];
        tdtxtfile = [tdsrcfld '/' tbase_opts.lfile];

        % initialize lookup volume
        lookup = uint8(0);
        lookup(tbase_opts.xsize, tbase_opts.ysize, tbase_opts.zsize, 5) = lookup(1);

        % read label file
        tbase_opts.labels = splittocell(asciiread(tdtxtfile), char([13,10]), 1, 1);

        % first line gives number
        labcount = str2double(tbase_opts.labels{1});
        tbase_opts.labels(1) = [];

        % remove emtpy lines at end (to make sure)
        while isempty(tbase_opts.labels{end})
            tbase_opts.labels(end) = [];
        end

        % checking file content
        if isempty(labcount) || ...
            length(tbase_opts.labels) ~= labcount

            % clear memory and bail out
            clear tbase_opts;
            error( ...
                'BVQXtools:IllegalFileContent',...
                'Illegal TD label file found.' ...
            );
        end

        % parse labels
        for lc = 1:labcount

            % split at comma and parse entry number
            tbase_opts.labels{lc}    = splittocell(tbase_opts.labels{lc}, ',', 1);
            tbase_opts.labels{lc}{2} = str2double(tbase_opts.labels{lc}{2});
        end

        % binary read source file
        try
            fid = fopen(tdsrcfile, 'r');
            tdlen = fread(fid, [1, 4], 'uint8=>double');
            tddbc = fread(fid, Inf, '*int8')';
            fclose(fid);
        catch
            clear tbase_opts;
            error( ...
                'BVQXtools:ErrorOpeningFile',...
                'Database file not found/unreadable: ''%s''.',...
                tdsrcfile ...
            );
        end

        % get number of bytes and points
        numb = length(tddbc);
        numpoints = ...
            tdlen(1) * 16777216 + ...
            tdlen(2) *    65536 + ...
            tdlen(3) *      256 + ...
            tdlen(4);

        % tell what we are doing
        disp('No TD cache file found, building internal DB...');
        progresscount(numpoints, 4, 20, '.');

        % start at byte pos after length of points
        rpos = 1;
        nump = 0;
        tx   = tbase_opts.xtrans;
        ty   = tbase_opts.ytrans;
        txyz = [tx, ty, 0];

        % loop until all bytes processed
        while rpos < numb

            % get size of element
            esize = double(tddbc(rpos));

            % get coords, z representation and label index list
            coord = double(tddbc(rpos+1:rpos+3)) + txyz;
            zreps = find(tbase_opts.zvals == coord(3));
            ilist = double(tddbc(rpos+4:rpos+esize));
            ilist(ilist < 0) = 127 - ilist(ilist < 0);
            rpos  = rpos + 1 + esize;
            nump  = nump + 1;

            % show progress
            if ~mod(nump, 1000)
                progresscount(nump);
            end

            % check if z is OK
            if isempty(zreps)
                warning( ...
                    'BVQXtools:BadFileContent',...
                    'Illegal coordinate, not Talairach: %d, %d, %d',...
                    coord(1), coord(2), coord(3) ...
                );
                continue;
            end

            % check label index list
            if any(ilist > labcount)
                warning( ...
                    'BVQXtools:BadFileContent',...
                    'Illegal label list for coordinate %d, %d, %d',...
                    coord(1), coord(2), coord(3) ...
                );
            end

            % fill label lookup
            coord(3) = zreps(1);
            lpos = 1;
            while length(ilist) > 1
                if ilist(2) ~= 1
                    warning( ...
                        'BVQXtools:BadFileContent',...
                        'Illegal label list entry in point %d.', ...
                        nump ...
                    );
                    break;
                end
                lookup(coord(1), coord(2), coord(3), lpos) = uint8(ilist(1));
                ilist(1:2) = [];
                lpos = lpos + 1;
            end
        end

        % check if everything went right
        if nump ~= numpoints
            clear tbase_opts;
            error( ...
                'BVQXtools:ErrorOpeningFile',...
                'Database contains wrong number of points.' ...
            );
        end

        % else populate lookup
        tbase_opts.lookup = lookup;
        clear lookup;

        % clear no longer needed memory
        clear tddbc;
    end

    % pure labellist
    labcount = length(tbase_opts.labels);
    plabels = cell(1, labcount);
    for lc = 1:labcount
        plabels{lc} = lower(tbase_opts.labels{lc}{1});
    end
    tbase_opts.plabels = plabels;

    % write file?
    if tbase_opts.reload && ...
        exist(tdmatfile, 'file') == 0
        disp(' ');
        disp('TD caching enabled, saving internal DB cache...');
        try
            saved        = struct;
            saved.labels = tbase_opts.labels;
            saved.lookup = tbase_opts.lookup;
            if mainver < 7
                save(tdmatfile, 'saved');
            else
                save(tdmatfile, 'saved', '-v6');
            end
            clear saved;
        catch
            warning( ...
                'BVQXtools:ErrorSavingMATfile',...
                'Couldn''t save intermediate file for database in folder ''%s''.',...
                strrep(tdsrcfld,'\','\\') ...
            );
        end
    end
end

% performing argument parsing
if ischar(rtype)
    rtype = str2double(rtype);
end
if nargin > 3 && ...
    ischar(xc) && ...
    ischar(yc) && ...
    ischar(zc)
    try
        xc = str2double(xc);
        yc = str2double(yc);
        zc = str2double(zc);
    catch
        error( ...
            'BVQXtools:BadArgument', ...
            'Invalid char argument given.' ...
        );
    end
end
if nargin > 4 && ...
    ischar(cs)
    cs = eval(['[' strrep(cs, '#', ',') ']'], '[7,3]');
end
if rtype < 0
    rtype = -rtype;
    btype = 1;
else
    btype = 0;
end

% default set label components (empty) and return string (*)
rstring  = '*';
rxlabelc = [];
xlabel   = 'l_empty';
xlabelc  = zeros(1, 5);

% different lookup types
switch (round(rtype))

    % Talairach label
    case {2}

        % get lookup x,y,z coordinates and sizes
        xc = round(xc(1) + tbase_opts.xtrans);
        yc = round(yc(1) + tbase_opts.ytrans);
        zc = floor(zc(1));
        xm = tbase_opts.xsize;
        ym = tbase_opts.ysize;
        zn = tbase_opts.zmin;
        zm = tbase_opts.zmax;

        % check coordinates
        if  xc < 1 || ...
            xc > xm || ...
            yc < 1 || ...
            yc > ym || ...
            zc < zn || ...
            zc > zm
            return;
        end

        % find good z representation slice
        zc = tbase_opts.zlui(zc + 1 - zn);

        % get label index list
        lilist = double(squeeze(tbase_opts.lookup(xc, yc, zc, :)));
        lilset = lilist(lilist ~= 0);

        % any content
        if ~isempty(lilset)

            % reset output
            rstringc = {'*','*','*','*','*'};
            rxlabelc = [0 0 0 0 0];
            for lnum = lilset(:)'
                rstringc{tbase_opts.labels{lnum}{2}} = tbase_opts.labels{lnum}{1};
                rxlabelc(tbase_opts.labels{lnum}{2}) = lnum;
            end
            if btype
                rstring = {rstringc{1:5}};
            else
                rstring = [ ...
                    rstringc{1}, ',', ...
                    rstringc{2}, ',', ...
                    rstringc{3}, ',', ...
                    rstringc{4}, ',', ...
                    rstringc{5}];
            end
        end

    % Talairach cube lookup
    case {3}

        % only valid with 5 arguments, cs being double
        if nargin < 5 || ...
           ~isa(cs, 'double') || ...
            numel(cs) ~= 1
            error( ...
                'BVQXtools:TooFewArguments',...
                'Argument ''size'' missing.' ...
            );
        end

        % get good cs from [3..13]
        cs  = min(13, max(3, cs));
        cs  = floor(cs / 2);

        % create data holding structs
        lss = struct;
        lsn = struct;

        % iterate over all coordinates
        for xr = (xc-cs):(xc+cs)
            for yr = (yc-cs):(yc+cs)
                for zr = (zc-cs):(zc+cs)

                    % recursively call
                    [lc, lcf] = tdlocal(2, xr, yr, zr);

                    % check struct
                    if ~isfield(lss,lcf)
                        lss.(lcf) = 1;
                        lsn.(lcf) = lc;
                    else
                        lss.(lcf) = lss.(lcf) + 1;
                    end
                end
            end
        end

        % get filled in, distinct labels, init to zero (to sort by number)
        lt  = fieldnames(lsn);
        lti = zeros(size(lt));

        % iterate over fields (unique labels)
        for lout = 1:length(lt)
            lti(lout) = lss.(lt{lout});
        end

        % sort and init return string
        [ltii{1:2}] = sort(lti);
        ltii = ltii{2};
        rstring = '';

        % fill return string
        for lout = length(lt):-1:1
            rstring = sprintf('%s%d:%s:', rstring, ...
                lss.(lt{ltii(lout)}), ...
                lsn.(lt{ltii(lout)}));
        end

    % Nearest gray matter search
    case {5}

        % only valid with 5 arguments, cs being double
        if nargin < 5 || ...
           ~isa(cs, 'double') || ...
            numel(cs) ~= 2 || ...
            any(isinf(cs) | isnan(cs) | cs < 1 | cs > 10);
            error( ...
                'BVQXtools:TooFewArguments',...
                'Argument size missing or invalid.' ...
            );
        end

        % get good range for sphere
        xcs      = min(8.501, max(1.19, cs(1)));
        radius   = fix(xcs + 0.8);
        middle   = radius + 1;
        diameter = 2 * radius + 1;
        found(1:diameter, 1:diameter, 1:diameter) = 0;

        % check direct coordinate first
        tresult = tdlocal(2, xc, yc, zc);
        rstring = '*:*';

        % if found
        if strfind(lower(tresult), 'gray matter')
            if btype
                tresult = tdlocal(-2,xc,yc,zc);
                rstring = {{[xc,yc,zc], [0,0,0], 0, tresult}};
            else
                rstring = ['+0,+0,+0  (dist=0.000): ' tresult];
            end
            return;
        end

        % iterate over positions
        for xr = -radius:radius
            for yr = -radius:radius
                for zr = -radius:radius

                    % check distance to center
                    dist = sqrt(xr*xr + yr*yr + zr*zr);
                    if dist > xcs
                        continue;
                    end

                    % call tdlocal recursively
                    tresult = tdlocal(2,xc + xr, yc + yr, zc + zr);

                    % if gray matter found
                    if strfind(lower(tresult), 'gray matter')
                        found(middle + xr, middle + yr, middle + zr) = dist;
                    end
                end
            end
        end

        % no coordinate found
        rfound = find(found(:));
        if isempty(rfound)
            return;
        end

        % find distances
        rdist  = found(rfound);
        [rsort,rindex] = sort(rdist);

        % how many distinct labels at most?
        cs(2) = fix(cs(2) + 0.1);

        % init return string
        if btype
            rstring = {};
        else
            rstring = '';
        end

        % init struct
        lss = struct;

        % iterate over number of found occurrences
        radius = radius + 1;
        for rcount = 1:length(rsort)

            % get coordinate for index
            [rpx, rpy, rpz] = ind2sub(size(found), rfound(rindex(rcount)));
            mxc  = xc - radius + rpx;
            myc  = yc - radius + rpy;
            mzc  = zc - radius + rpz;

            % on requested output type
            if btype

                % call tdlocal again
                [tresult, lcf] = tdlocal(-2, mxc, myc, mzc);

                % put into output
                if ~isfield(lss, lcf)

                    % add to cell output
                    rstring{end+1} = {[mxc, myc, mzc], ...
                                      [rpx - radius, rpy - radius, rpz - radius], ...
                                       rsort(rcount), tresult};

                    % make field available and reduce cs(2)
                    lss.(lcf) = 1;
                    cs(2) = cs(2)-1;
                end

            % direct output
            else

                % call tdlocal again
                [tresult, lcf] = tdlocal(2, mxc, myc, mzc);

                % put into output
                if ~isfield(lss, lcf)
                    rstring = sprintf('%s%+d,%+d,%+d  (dist=%.3f): %s (coords=%d,%d,%d):', ...
                                      rstring, rpx - radius, rpy - radius, rpz - radius, ...
                                      rsort(rcount), tresult, ...
                                      mxc, myc, mzc);

                    % make field available and reduce cs(2)
                    lss.(lcf) = 1;
                    cs(2) = cs(2)-1;
                end
            end

            % leave if number of hits was found
            if cs(2) < 1
                break;
            end
        end

    % Talairach label search
    case {6}

        % prepare output
        rstring = [zeros(0,1), zeros(0,1), zeros(0,1)];

        % check second arg
        if ~ischar(xc) || ...
            isempty(xc)
            return;
        end

        % find matching strings
        ccs = lower(xc);
        ccm = find(strcmp(tbase_opts.plabels, ccs(:)'));
        if isempty(ccm)
            return;
        end

        % find matching voxels
        mvx = [];
        for lc = ccm
            mvx = union(mvx, find(any(tbase_opts.lookup == uint8(lc), 4)));
        end
        [mvx, mvy, mvz] = ind2sub( ...
            [tbase_opts.xsize, tbase_opts.ysize, tbase_opts.zsize], mvx);

        % create output array
        mvx = mvx - tbase_opts.xtrans;
        mvy = mvy - tbase_opts.ytrans;
        ccm = sort(unique(mvz(:)));
        tx = zeros(0, 1);
        ty = zeros(0, 1);
        tz = zeros(0, 1);
        for lc = ccm(:)'
            zmi = find(mvz == lc);
            zmc = tbase_opts.zlub(tbase_opts.zlui == lc);
            for zc = zmc(:)'
                tx = vertcat(tx, mvx(zmi));
                ty = vertcat(ty, mvy(zmi));
                tz = vertcat(tz, zc * ones(length(zmi), 1));
            end
        end
        rstring = [tx,ty,tz];

    % invalid type
    otherwise
        rstring = '* (unknown rtype in request!)';
end

% add second output argument if requested
if nargout > 1
    if ~isempty(rxlabelc)
        xlabel = sprintf('l_%d_%d_%d_%d_%d',rxlabelc);
    end
    if nargout > 2
        if ~isempty(rxlabelc)
            xlabelc = rxlabelc;
        end
    end
end
